A simple and effective Verlet-type algorithm for simulating Langevin dynamics 
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We present a revision to the well known Stormer-Verlet algorithm for simulating second order dif- 
ferential equations. The revision addresses the inclusion of linear friction with associated stochastic 
noise, and we analytically demonstrate that the new algorithm correctly reproduces diffusive be- 
havior of a particle in a flat potential. For a harmonic oscillator, our algorithm provides the exact 
Boltzmann distribution for any value of damping, frequency, and time step for both underdamped 
and overdamped behavior within the usual stability limit of the Verlet algorithm. Given the struc- 
ture and simplicity of the method we conclude that this approach can trivially be adapted for 
contemporary applications, including molecular dynamics with extensions such as molecular con- 
straints. 

PACS numbers: 



I. INTRODUCTION 

In Molecular Dynamics (MD) simulations, Newton's 
equations of motion are solved numerically to produce 
trajectories of the system in the microcanonical [N, V, E) 
ensemble. The most commonly used scheme for this pur- 
pose is the Verlet method, which is based on truncated 
Taylor expansions for the evolution of a particle with 
mass m, coordinate r(i), velocity and force f{r,t) 
Introducing the discrete-time variables r" = r(t„), 
= v{tn), and /" = f{r",tn), the so-called velocity ex- 
plicit Verlet (or velocity Verlet) scheme reads 
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By considering two successive time steps and eliminating 
the velocity variables from the equations, the more origi- 
nal form of the Stormer-Verlet method [1, 0] (here called 
the position- Verlet method), is found 
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with the associated velocity calculated by the central dif- 
ference 
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The Verlet scheme is accurate to second order in dt, which 
means that the deviation, per time step, of the computed 
trajectory from the true (analytic) one, scales with the 
third power of dt. The merits of the Verlet scheme that 
make it so widely popular for MD simulations include 
its convenience, efficiency, and time reversibility, which 
ensures that the error of the total energy of long-time 



integrations is bounded and does not drift or diffuse, as 
illustrated by the exact solution to the Verlet method 
applied to a harmonic oscillator (see, e.g., below). 

The most frequently used ensemble in statistical- 
mechanics is the canonical {N, V, T) ensemble where the 
temperature of the system, rather than its energy, is con- 
stant. A variety of methods for conducting MD simula- 
tions in the canonical ensemble have been proposed over 
the years. A very appealing class of such methods include 
integrators for Langevin Dynamics (LD) simulations. In 
LD, two forces are added to the conservative force field - 
a friction force proportional to the velocity with friction 
coefficient a > 0, and thermal white noise /3(i). Elxplic- 
itly, the Langevin equation of motion is given by j5| 



r = V 
mv — f{r, t) ^ av + j3{t) . 
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In order to satisfy the dissipation-fluctuation theorem, 
it is often assumed that the stochastic force is Gaussian 
distributed, and has the statistical properties Q 
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where fcs is Boltzmann's constant and T is the thermo- 
dynamic temperature of the heat bath. 

The difficulty to develop accurate integrators for LD 
stems from the non-analytic nature of /3(i), which in- 
validates the Taylor expansion commonly used for the 
derivation of the Verlet scheme. The most naive way to 
overcome this difficulty is to replace the delta-function 
correlated noise with a set of rectangular pulses of mean- 
squared size y^2kBTa/dt, each of which acting over the 
centered time interval (t„ — dt/2,tn + dt/2). Employing 
this approximation for (3{t) yields the famous Brooks- 
Briinger-Karplus (BBK) scheme 0, which unfortunately 
turns out to exhibit a simulated temperature that differs 
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by 0{dt) from the correct one 0- This disappointing 
result can be attributed to the effect of a small time step 
dt, which results in a stochastic velocity change of or- 
der \/di that overshadows the "regular" (deterministic) 
linear velocity change of order dt. Statistical analysis of 
the small time interval shows, in fact, that the opposite 
is true, since the average over all noise realizations van- 
ishes. In order to develop a reliable integrator for LD 
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one therefore needs to carefully treat the coupling be- 
tween the stochastic and analytic contributions. This is, 
unfortunately, not easily accomplished with a discretized 
approximation for l3{t). A different approach to the prob- 
lem has been introduced by Ermak and Buckholz (EB) 
[8j . The EB method is based on a numerical integration 
of the formal solution to the Langevin equation. This 
gives 



= exp(^adt/m) + — f exp [-a(i„+i - r)/m] [f{t') + /3(i')] dt' (9) 

rnv" 1 
r»+i = [1 - exp {-adt/m)] + - {1 - exp [-a(i„+i - t')/m]} [f{t') + p{t')] dt' . (10) 



Several features distinguish the EB method from Verlet- 
type schemes such as BBK. First, the friction force is 
accounted for separately from the other two forces, via 
exponentially decaying functions. Second, this scheme 
considers an ensemble of trajectories instead of a single 
one, thereby requiring two random Gaussian numbers per 
time step. The one appearing in the velocity equation, 
/jj*^* exp[— Q;(c?t — t')/m]l3{t')dt' , represents the stochastic 
change in velocity during the time interval. The other, 
— exp[—a{dt — t')/m]}l3{t')dt' , appears in the po- 
sition equation and characterizes the random displace- 
ment. These two numbers are different, but correlated. 
Third, there is formally no need to compute the new po- 
sition before the new velocity as in the velocity- Verlet 
scheme and, in principle, one can interchange the order 
by which the equations are evolved. The last point is 
particularly important since it is directly related to the 
most critical issue that remains unsolved with the EB 
approach, which is how to perform the numerical inte- 
gration over the deterministic force f{t). The appear- 
ance of exponential "weight functions" in the integrals in 
Eqs. © and ([TU]) opens the possibility for using a vari- 
ety of linear combinations involving /"~^, /", and /"''"^. 
The most popular integration schemes for the determinis- 
tic force constitute the van Gunsteren-Berendsen Q and 
the Langevin impulse [l^l methods. 

More recently, Ricci and Ciccotti (RC) introduced a 
formalism for a systematic derivation of numerical inte- 
grators for LD The essential steps of the formalism 
are to (i) express the integrator as an exponential oper- 
ator, (ii) split the time interval dt into a series of even 
smaller time steps, and (iii) use the Suzuki- Trotter ex- 
pansion for time-ordered exponential operators to find 
an approximation for the evolution operator correspond- 
ing to dt. By using mid-point splitting, RC arrived at 
a scheme that require only one random number per dt. 
The same formalism was later used to develop a new fam- 
ily of integrators that more carefully treat the coupling 



between the stochastic and deterministic components of 
the dynamics [l^, [Si • ^ characteristic feature of these 
new schemes is that they require two independent ran- 
dom numbers per time step, which is different from the 
EB family of schemes, where two correlated random num- 
bers are involved. 

Returning for a moment to MD simulations of Newto- 
nian Dynamics, it is important to recall that, due to its 
computational efficiency and time reversibility, the Verlet 
method is regarded as superior to other integration meth- 
ods that may produce more accurate trajectories. In the 
case of LD, the trajectories include a random component 
and, therefore, their accuracy (which can only be defined 
in statistical terms) becomes an issue of even less impor- 
tance. The efficiency of different methods for canonical 
ensemble simulations of many-particles systems must be 
tested according to their ability to reproduce measurable 
statistical quantities, such as the Boltzmann distribution 
corresponding to the temperature of the heat bath. To 
this end, none of the existing methods has managed to 
demonstrate exact thermodynamic response ^Ti] . We here 
present a novel Verlet-type scheme for LD simulations 
that is very simple to implement and which yields cor- 
rect statistical-mechanical behavior of a particle diffusing 
in both a flat and a harmonic potential. 



II. DERIVATION OF THE ALGORITHM 



In the spirit of the simplicity of the Verlet algorithm 
in its traditional forms, we here arrive at a new useful 
method through a straight-forward derivation. Starting 
with the continuous-time Langevin equations ([S])-®, we 
integrate Eq. ^ over a (small) time interval dt between 
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two times, t„ and i„+i = in + dt: 



mv dt' ^ f dt' ~ ar dt' 
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P{t')dt' . (11) 



With no approximation, this can be written 



An equation similar to (ITSI) has been used by Ricci and 
Ciccotti in one of the forms of their scheme (see Eq.(18) 
in Ref. [HI). It introduces errors that scale with dt'^ in 
the deterministic trajectory and in the variance of the 
stochastic component (see discussion in Ref. 14]). 

Inserting ?;"+^ from Eq. (fT2|) into Eq. ([TSl) provides a 
convenient pair of equations 



m(t)"+i -«")=/ / dt' - a(r"+i - r") + /3" 



where 
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fi{t') dt' 
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is a Gaussian random number with (/3") = and 
(^"^') = 2akBTdtSnj. 
Integrating over Eq. (O yields 



rdt' = r"+i -r" = 

which can be approximated with 
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where 



b = 



1 I grf^ 
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For any given realization of /3", Eq. ([T5| is correct to 
second order in dt, while (fT7| is exact as written. We 
now approximate the integral of the deterministic force 
/ such that both equations are correct to second order in 
dt: 
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We first notice that when a = 0, the above equations 
(IT9l) - (|20l) reduce to the standard velocity- Verlet scheme 
given in Eqs. ([1]) and ([2]). Second, we notice the very 
reasonable feature of the representation of the damp- 
ing, which is calculated as the integral of the actual 
path that the object has traveled during the time step 
dt. Third, the noise is represented as a single stochas- 
tic variable for each time step, realized by a single ag- 
gregated impulse during dt that influences the dynam- 
ics over the time step. In this regard, Eqs. (IT9|) and 
(PO)) constitute a simple functional Verlet-type scheme for 
solving stochastic Langevin equations. Unlike the afore- 
mentioned EB-type family of schemes (8l-[l0| , our method 
does not employ two stochastic variables. This is a con- 
sequence of our main objective; namely to produce the 
correct statistical-mechanics for large ensembles (spatial 
or temporal) , while focusing less on the detailed dynam- 
ics within each time step. 



Before proceeding to analyzing the behavior of the de- 
veloped scheme, we rewrite the method in a couple of 
equivalent different forms that can be useful and that 
can illustrate close connections to previously published 
work. First, we observe that by inserting Eq. into 
([20|. the equations can be rewritten 



r"+i = r'- + bdtv^ + —r + —r+^ 
2m 2m 
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This form reveals that the derived method presented 
here parallels one mentioned in passing by Melchionna 
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(Eq. (41) in Ref. where the contribution of the fric- 

tion force has been integrated similarly, but with a dif- 
ferent noise term that includes two independent random 
numbers per time step. Second, we rewrite the scheme in 
the form given only by displacement coordinates. This 
is accomplished by subsequently inserting the expressions 
for w" ^ and r" ^ into that of r"+i given in Eq. (CT . 
The result is: 



m 2m 
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with the associated velocity given by 
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Eq. (|24p shows that our formulation is also in close 
proximity to the BBK scheme Q , which is of a position- 
Verlet-type. (Notice that it is essential to have the proper 
form for the accompanying velocity in order to determine 
if a position formulation is the same as a velocity explicit 
form.) However, unlike our Eq. the BBK scheme 

employs only a single stochastic number for the two time 
steps covered in the displacement equation. In fact, the 
methods of Melchionna and BBK, while not identical, are 
closely related, since the use of a single stochastic num- 
ber in the position representation (BBK) translates into 
applying two random numbers in the velocity equation 
(Melchionna). Conversely, in our scheme, we have a sin- 
gle random number in the velocity formulation, and two 
in the position representation covering two time steps. 



III. LINEAR ANALYSIS 

In order to evaluate the general applicability of the 
above method, we calculate key statistical measures of 
the numerical scheme for some characteristic linear cases 
p^ . and compare them to known statistical- mechanical 
behavior of the true Langevin system. 



A. Thermal Diffusion, / = 

The first basic property to investigate is the diffusive 
behavior of a particle moving in a fiat potential (/ — 0) 
at temperature T. In this case, Eq. reads: 
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By using the same equation to express in terms of v"^ ^ 
and and by repeating this procedure until reaching 



the initial velocity u", one arrives at the relation 
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where cr is a standard Gaussian random number with 
(ct) = and (cr^) — 1, the superscript denoting differ- 
ent realizations of this variable such that {a" a'') — 5ni- 
Summing over the random numbers in Eq. (j28p yields 
another Gaussian random number, and in combination 
with Eq. ([T5)) . one arrives at 
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For large n, beyond the transients from the initial condi- 
tion (a" <C 1), we find that the velocity is characterized 
by a Gaussian (Maxwell-Boltzmann) distribution with 
zero mean and 
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which results in reproducing the exact expectation for 
the average kinetic energy (thermal energy) 
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Complementing this result, we turn to the displacement 
coordinate Eq. for / = 0: 
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Inserting Eq. (^5]) for ^ , and repeating the procedure 
until reaching the initial position r°, we obtain the result 
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By using the defintion of a and b ((T5)) . we find that 
for large n (such that a" <^ 1), Eq. can be written 
as 
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The second term in Eq. f34|) . which is the transient bal- 
listic displacement, matches exactly the corresponding 
value predicted by the Langevin solution [second term 
on the r.h.s of Eq. (ITOl) for — >■ 00]. More importantly, 
we find that our scheme results in a simulated diffusion 
coefficient 
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which agrees with Einstein's fluctuation-dissipation re- 
lationship for LD j5i]. Notice that the exact results for 
fluctuations and diffusion obtained here are independent 
of the magnitudes of the time step dt, damping a > 0, 
and temperature T. 



B. Thermal Harmonic Oscillator, / — —Kr 

Encouraged by the performance of the method for 
/ = 0, we now turn to analyzing the method applied to 
a damped thermal harmonic oscillator, where f — —nr 
represents a linear Hooke's spring with spring constant 
K > 0. Our Eqs. ^ and dH]) now read 
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and where Q,q = ^ K/m is the resonance frequency of the 
undamped continuous-time oscillator. The matrix V has 
the eigenvalues K±: 
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(39) 



from where we can evaluate the formal stability limit 
VLi^dt < 2, consistent with the requirement |A±| < 1. 
Notice that A_|_A_ = a for all parameter values. 

We analyze the basic thermodynamic properties of the 
thermal harmonic oscillator by perpetuating the map 
Eq. (|36p from initial conditions r*^ and v'^, 



= V 
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regimes of simulated oscillator behavior as a function of 
those two parameters. The two physically relevant dy- 
namical regimes, underdamped (A) and overdamped (B), 
are considered first. 
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FIG. 1: Sketch of the characteristic regimes for the Verlet in- 
tegrator applied to the damped harmonic oscillator. Regimes 
are defined by the eigenvalues given in Eq. (|39|l . Solid curve is 
the boundary between physical (left) and unphysical (right) 
behavior. Vertical dash-dotted line is the formal stability 
boundary of the Verlet method, A„ = —1. Marker • in- 
dicates where A± = 0. Regimes are: (A) Underdamped, 
A± complex; (B) Overdamped, < a < A± < f real; (C) 
— f < A± < —a < real; (D) A± real, A_ negative, A+ 
positive (o < 0); (E) Unstable, A_ < -f . 



me A, Underdamped Dynamics 



The oscillator is always stable in this regime, where 
(a/2m)^ < ilo(l — f^o'^i^/4), and the eigenvalues can be 
conveniently written 



Ah 



(43) 



k=0 



(41) Here, fiy is the underdamped Verlet-oscillator resonance 
frequency given by 



where the fcth iteration is conveniently given by the uni- 
tary transformation 
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The two characteristic parameters of the Verlet 
method applied to a damped harmonic oscillator are 
a/2mQ.Q and ^^dt. Figure [T] displays the five different 
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■y/a cos riydt + 1 — 6 
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\fa cos fiy + a — & 
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The diagonalizing transformation is then given by 



\u\ 



+ sinfiydi i^sinfiydt 



(47) 



with |m| being the normahzation of the eigenvectors of V that appear as the column vectors in U. The fcth time step 
can then be expressed from Eq. (|^^ in the exphcit form: 



bdt 



^/a sin fly dt 



n-k ( ^ singly (ii COS riyfcdt + ^ sin ilyfcdi sinily/cdi \ 

^ I -(p^sin^rjydt + (^)2)sinr2yfcdt ^ sin fiydt cosily fcdt - 3^ sin rjyfcdi ] 



from where we obtain the explicit form of the expression of interest in Eq. (j4T 
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) -^[y/a sin fly dt cosily kdt — (1 — ^/a cos ^ly dt) sinilykdt] 
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The thermodynamic limit of ((r")^) and ((w")'^) can now summation is an independent Gaussian random number, 
be found by inserting Eq. (|49| into Eq. (j4T|) . whereafter Taking the limit n — cx) this yields for the variance of 
the square sum can be evaluated since each term in the the displacement: 



J 



((r")2) = 



aksTb'^dt^ 
2m?a siv? ilydt 



-(1 + a + cos fly dt) ^ ( 



fe=0 
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-(1 + a cos 2nydt + 2-v/acos flydt) a'^ cos 2ilykdt 



k=0 



sin rtydt{l + cos ftydt) a'^ sin 2VLykdt 



k=0 



aksTb'^dt^ 
2m?'a siv? ilydt 



146(1-^^) (1 - a)Vacosrjydt + i(l - a^) 
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1 - a 
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We similarly obtain the result for the statistical variance of the velocity: 
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2akBTlPdt^ 
TO^asin^ flvdt 
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1 + — 2a cos 2il,vdt 

'odt^ \ _ ( a \2 
4 ' y2m) 



1 - 
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The two variances are noteworthy, since evaluating the 
average potential and kinetic energies gives 

= \^{{rn^) = \kBT (52) 
Ek = \m{{v- f) = hsT [l ^) . (53) 

These hold true for any parameter choice in the un- 
derdamped regime. Notably, Eq. ([52)) implies that our 
method produces the exact statistical distribution re- 
gardless of time step, frequency (potential curvature), 
damping parameter, or temperature. Therefore, it is 
reasonable to expect that this method will provide cor- 
rect Boltzmann distribution in thermodynamic equilib- 
rium when simulating complex many-particle systems 
that may have a multiple of participating frequencies. To 
our knowledge, this important feature has not previously 
been reported for a numerical integrator of the Langevin 
equation [7]. 

Complementing the Boltzmann distribution of the dis- 
placement is the variance of the velocity, Eq. (|53)) . De- 
spite the obvious discrepancy between the true kinetic 
energy of the Langevin equation and the one shown for 
our algorithm, we submit that the presented result is the 
best possible for a method that builds on the discretized 
Verlet formalism. The reduction in the variance of the 
velocity by a factor of (1 — Vi^dt^ /A) does not arise from 
the treatment of friction and noise, as implied by the fact 
that this factor depends on neither damping nor temper- 
ature. Instead, the observed deviation arises from the 
approximation of the potential curvature introduced by 
the discrete integrator. This is consistent with the well- 
known inherent discrepancy between displacement and 



associated velocity that causes the periodic deviations 
with magnitude dt^ from strict energy conservation in 
a simulated harmonic oscillator. This can be explicitly 
demonstrated by using the Verlet equations (H]) and ^ 
for an undamped harmonic oscillator with initial condi- 
tions r'^ and The result is 




(54) 

which illustrates that the Verlet velocity (momentum) is 
depressed and is not exactly the conjugated coordinate 
to the displacement, and that the discrepancy is related 
to the proportionality seen between the two thermody- 
namic expressions ([5^ and (15^ . 



2. Regime B, Overdamped Dynamics 

This regime is defined by the requirements (a/2m)'^ > 
- nidi's/A) and adt/2m < 1 ( for a > 0). The 
latter condition is imposed to ensure A± > 0, which is 
necessary for physically meaningful dynamics, where only 
monotonic decay is possible in the overdamped regime for 
T = 0. 

The matrix V appearing in Eq. p6l) . which is now ex- 
pressed by 
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has the eigenvalues 
where Ay is given by 



, ^/a cosh Xv dt + 1 — b bdt \ 

I ^Mf[(^)2 _ _^sinh2 Aydt] y/^coshXvdt + a-b J ' 



A± = V^e*^^''*, (56) 



V^coshXvdt = b(l-^^] (57) 



V^sinh Aydt = bdt J - f^o ( 1 ^ ) ■ (58) 



The matrix V can be diagonahzed (see Eq. p6p) using the transformation matrix 

" ^ ( +"^sinhAWt] ^[-^ J"isinhAydt] ) ^^^^ 

with and being the normahzations of the two real eigenvectors of V. This leads to the explicit results: 

^f. _ bdt r-k I ^ sinh Xydt cosh Xykdt + ^ sinh Xykdt sinhAyfcdf \ 

V^sinhAydt " _ (^(^^f _ si^^i^"^ Xy ^t) sinh Xy kdt smh Xydt cosh Xy kdt - ^sinhXy kdt J 

(60) 

and 



^,kKron~k bdty/2akBT _fe 

vW" = - — ^ . , ■ X 

Zmy/a smh Aydr 
Y^sinh Aydt coshAyfcdt + (1 + ^/acosh Aydi) sinhAyfcdt 
_ ^lo^'j 2_^^g:^^j^xYdt cosh Xykdt — (1 — yacoshAydt) sinhAyfcdt] 

I 



(7"-'=. (61) 



Following the same procedure used in the underdamped 
case above, we insert Eq. (pTj) into Eq. (PT|) . and eval- 
uate the variances of the displacement and velocity in 
the thermodynamic limit for n — )■ oo by calculating the 



square sum of the amplitudes of the independent stochas- 
tic numbers of each term. After some calculations, we 
arrive at 



akeTb^dt^ [ 1., o ^ i x ^ fe 

- -(1 + a + 2VacoshAydt)2_^a'' 



1rn?a s\vA\ Xydt 



fc=0 



^ CO 

-(1 + acosh2Ay(it + 2-v/acosh Aydt) a'^ cosh2Xvkdt 

k=0 

oo 

sinh Aydt (1 + 2A/a cosh Xydt) a*" sinh 2Xvkdt 



k=0 



aksTb'^dt^ 
2'm'^a sinh^ Ayrfi 

akeTb'^dt'^ 
2m?a sinh^ Xydt 
ksT 



146(1-^^) (l-a)VacoshAydi+i(l-a2)- 



2 1 -a 

2m "O 



1 -I- a2 — 2acosh2Aydi 



(62) 
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which is exactly the same result as for the underdamped case Eq. (|50p . We similarly obtain the result for the variance 
of the velocity: 



2akBTb'^dt 
T?a sinh'^ Xydt 



— -^(1 + a — 2Vacosh Ayd<) ^ i 



fc=0 



^ CO 

--(1 + acosh2Av'(it — 2-y/acosh Aydt) a*^ cosh2Ayfc(it 

fe=0 

oo 

-\/asinh A\/(i<(l — -yacosh Xydt) a'^ smh2Xv kdt 



2akBTb'^dt 
rn^a sinh^ Xydt 

2akBTb'^dt^ 
rn^a sinh^ Xydt 



2^+2 \ 2 ^^ illdt^ 



n^dt 



(1 — a)-ya cosh Ay (it — i(l — 



1 -a 



1 - 



1 + — 2acosh2Ayfii 

4 / V 2m / 



4— 

^2m 



(63) 



which is also identical to the comparable underdamped 
result from Eq. (|5ip . The averages of potential and ki- 
netic energies are therefore given by Eqs. ((5^ and ([55)) 
also for the overdamped regime. Thus, we can now con- 
clude that our method produces the exact Boltzmann 
distribution for any physical dynamics, underdamped or 
overdamped, regardless of time step, frequency (potential 
curvature), damping parameter, or temperature. 

3. Regime C, Q.v = Tr/dt 

This somewhat unphysical regime, is characterized by 
conditions similar to the overdamped regime B; i.e., 
adt/2m < 1 (for a > 0) and (^)2 > ^1(1-^^). How- 
ever, while regime B corresponds to il^dt < \/2, regime 
C is defined for flodt > v^, resulting in A± < 0. Thus, 
this regime is typically reached for large time steps when 
simulating lightly damped dynamics, and it is the pre- 
curser for violating formal stability, given by flgdi < 2 
(see Fig.[T|). At T = and with initial conditions r° 7^ 0, 
the dynamics is characterized by r" and alternating 
their signs every time step (i.e., fly — ir/dt) with an 
exponentially decaying envelope. For T > 0, the ther- 
modynamic properties of this regime can, therefore, be 
evaluated similarly to that of the overdamped regime, 
since we here can write 

= {-lfU^^+ I, '^U-\ (64) 

where U and A± are given by Eqs. (l5Bl) - ([5^ . Apart from 
the alternating sign change in Eq. ([64)) this is identical 
to the corresponding mapping in regime B. This implies 
that Eqs. ([5^ and ([55)) from regime B for the variances 
of r" and w" also apply here. Thus, we conclude that the 



thermodynamic results Eq. (15^ and Eq. ([55)) are found 
also in this regime, which is physically unreasonable, yet 
numerically accessible. 



4. Regime D, Unphysical 

This unphysical regime, adt/2m > 1 (yet within the 
formal stability criterion n^dt < 2), is typically reached 
for large time steps when simulating moderate to strongly 
damped dynamics (see Fig.jT]). Since this regime involves 
one positive and one negative eigenvalue, we can map 
this onto the analysis from regime B with ^/a y^joj 
and Xy — > —Ay to complete the analysis and obtain the 
Boltzmann distribution. 



IV. DISCUSSION AND CONCLUSION 

We have presented a new Verlet-type algorithm for 
simulating Langevin dynamics. The method, written 
in the three different forms Eqs. ([19]) -([20]), Eqs. ([2T))- 
([^ or Eqs. ([M])-((25]), is aligned with other published 
methods that have demonstrated first, second, and third 
order accuracy for both simulated trajectories and de- 
rived thermodynamic quantities [B-EH. Linear analysis 
demonstrates that the method of this paper is robust and 
capable of providing exact representation of both diffu- 
sive behavior in a flat potential and the Boltzmann distri- 
bution in a harmonic potential regardless of damping and 
frequency (curvature of the potential) . The derived exact 
distribution is obtained for any time step subject to the 
usual Verlet stability limit and the condition adt < 2m, 
which is a necessary requirement for a meaningful atten- 
uated trajectory. The method is very simple and in the 
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usual Verlet structure, which means that it can be read- 
ily implemented for any Langevin application, including 
molecular dynamics of many-particle systems with and 
without molecular constraints and other commonly used 
modeling features. 

As a final note, we underline that a Verlet-simulated 
oscillator may not measure the exact temperature from 
the variance of the velocity (see Eq. (IS^I) ). This de- 
ficiency, which is common to all Verlet-type methods, 
arises from the known discrepancy between displacement 
and momentum as conjugated variables. This means that 
using the variance of the velocity for precisely assessing 



the temperature of a simulated system may be counter- 
productive. The appropriate criterion for obtaining the 
desired temperature is the achievement of correct statis- 
tical sampling, which is indeed given by our method. 
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